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Disclaimer 


No  warranties 3  express  or  implied,  are  made 
by  the  distributors  or  developers  that 
STARPAC  or  its  constituent  parts  are  free  of 
error.  "Vh.ey  should  not  be  relied  upon  as 
the  sole  basis  for  solving  a  problem  whose 
incorrect  solution  could  result  in  injury  to 
person  or  property .  If  the  programs  are 
employed  in  such  a  manner,  it  is  at  the 
user's  ovm  risk  and  the  distributors  and 
developers  disclaim  all  liability  for  such 
misuse. 

Computers  have  been  identified  in  this  paper 
in  order  to  adequately  specify  the  sample 
programs  and         test         results.  Such 

identification  does  not  imply  recommendation 
or  endorsement  by  the  National  Bureau  of 
Standards,  nor  does  it  imply  that  the 
equipment  identified  is  necessarily  the  best 
available  for  the  purpose. 
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Preface 


STARPAC,  the  Standards  Time  Series  and  Regression  Package,  is  a  library 
of  Fortran  subroutines  for  statistical  data  analysis  developed  by  the 
Statistical  Engineering  Division  (SED)  of  the  National  Bureau  of  Standards 
(NBS),  Boulder,  Colorado.  Earlier  versions  of  this  library  were  distributed 
by  the  SED  under  the  name  STATLIB  [Tryon  and  Donaldson,  1978] .  STARPAC 
incorporates  many  changes  to  STATLIB,  including  additional  statistical 
techniques,  improved  algorithms  and  enhanced  portability. 

We  are  distributing  STARPAC  source  code  in  segments,  and  the 
documentation  as  a  series  of  NBS  Technical  Notes.  This  will  allow  us  to  make 
the  new  capabilities  of  STARPAC  available  more  quickly  and  to  facilitate  the 
release  of  future  changes. 

The  first  of  the  series  of  Technical  Note  documenting  STARPAC  is  an 
Introduction  to  STARPAC  [Donaldson  and  Tryon,  1983]  which  gives  an  an  overview 
of  the  STARPAC  library,  defines  conventions  used  to  specify  the  documentation, 
and  presents  general  background  material.  The  Introduction  to  STARPAC  incudes 
information  which  is  essential  for  using  the  STARPAC  library,  and  users  should 
be  familiar  with  its  contents  before  attempting  to  use  any  STARPAC 
subroutine. 

This  Technical  Note  documents  the  STARPAC  nonlinear  least  squares 
regression  subroutines.   These  subroutines: 

•  use  the  quasi-Newton  NL2S0L  algorithm  [Dennis  et  al.  ,  1981a] ,  designed 
specifically  for  problems  which  may  have  highly  nonlinear  models,  large 
residuals,  or  poor  scaling; 

•  accept  user-supplied  derivatives  (the  Jacobian  matrix)  when  available, 
and  approximate  the  derivatives  numerically  otherwise; 

•  provide  utility  procedures  to  verify  the  correctness  of  the 
user-supplied  derivatives,  or  to  select  optimum  step  sizes  for 
numerically  approximating  the  derivative  [Schnabel,  1982]; 

•  allow  both  weighted  and  unweighted  analyses; 

•  permit  subsets  of  the  model  parameters  to  be  treated  as  constants  with 
their  values  held  fixed  at  their  input  values,  allowing  the  user  to 
examine  the  results  obtained  by  estimating  subsets  of  the  parameters  of 
a  general  model  without  rewriting  the  model  subroutine; 

•  are  easy  to  use,  providing  three  levels  of  user-control  of  the 
computations  and  results,  extensive  error  handling  facilities, 
comprehensive  printed  reports,  and  no  size  restrictions  other  than 
effective  machine  size; 


•  are  portable;  and 

•  are  easily  used  with  other  Fortran  subroutine  libraries. 

Other  code  segments  of  STARPAC  include  subroutines  for  time  series 
analysis  (in  both  time  and  frequency  domains),  line  printer  graphics,  basic 
statistical  analysis,  linear  least  squares  regression,  and  nonlinear 
optimization.  Users  may  obtain  a  list  of  the  available  code  segments  and 
their  associated  documentation  from  their  local  installer  of  STARPAC,  or 
directly  from  the  author.  Notation,  format,  and  naming  conventions  will 
remain  constant  throughout  the  series  of  Technical  Notes  documenting  STARPAC, 
allowing  the  documentation  for  each  code  segment  to  be  used  alone  or  in 
conjunction  with  the  documentation  for  another. 

The  code  for  STARPAC  was  developed  at  the  U.S.  Department  of  Commerce 
Boulder  Laboratories  on  a  CDC  CYBER  170/750  under  NOS  1.4,  and  all  examples 
presented  in  this  documentation  were  executed  in  this  environment  using  the 
FTN  4.8+552  compiler  with  rounding.  STARPAC  has  also  been  tested  on  the 
following  equipment. 


Computer 


Facility 


IBM  4  340 
VAX  11/780 


National  Center  for  Atmospheric  Research, 
Boulder,  Colorado 


Sperry  1100/82 


National  Bureau  of  Standards, 
Gaithersburg,  Maryland 


Perkin-Elmer  3230 


National  Bureau  of  Standards, 
Boulder,  Colorado 


Data  General  ECLIPSE  S/140 


Forest  Fire  Laboratory, 
Riverside,  California 


STARPAC  is  written  in  Fortran  '66.  Adherence  to  a  portable  subset  of 
ANSI  Fortran  [1966]  has  been  verified  by  the  PFORT  Verifier  [Rvder,  1974]. 
Workspace  and  machine-dependent  constants  are  supplied  using  subroutines  based 
on  the  Bell  Laboratories  "Framework  for  a  Portable  Library"  [Fox  et  al. , 
1978a].  We  have  also  used  subroutines  from  UNPACK  [Dongarra  et  al.  ,  1979], 
from  the  "Basic  Linear  Algebra  Subprograms  for  Fortran  Usage"  [Lawson  et  al. , 
1979],  and  from  DATAPAC  [Filliben,  1977].  The  subroutines  used  to  compute  the 
least  squares  solution  are  those  referenced  in  Dennis  et  al.  [1981a  and 
1981b],  and  the  algorithms  used  to  select  optimum  step  sizes  for  numerical 
derivatives,  and  to  check  analytic  derivatives  were  developed  by  Schnabel 
[1982].  The  printed  reports  from  the  nonlinear  least  squares  subroutines  have 
been  adapted  from  OMNITAB  II  [Hogben  et  al. ,  1971]. 
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Nonlinear  Least  Squares  Regression 
Using  STARPAC 
The  Standards  Time  Series  and  Regression  Package 

Janet  R.  Donaldson  and  Peter  V.  Tryon 

National  Bureau  of  Standards 
Washington,  DC  20234 

STARPAC,  the  Standards  Time  Series  and  Regression  Package,  is  a 
library  of  Fortran  subroutines  for  statistical  data  analysis 
developed  by  the  Statistical  Engineering  Division  (SED)  of  the 
National  Bureau  of  Standards  (NBS),  Boulder,  Colorado.  Earlier 
versions  of  this  library  were  distributed  by  the  SED  under  the  name 
STATLIB  [Tryon  and  Donaldson,  1978].  STARPAC  incorporates  many 
changes  to  STATLIB,  including  additional  statistical  techniques, 
improved  algorithms  and  enhanced  portability. 

STARPAC  emphasizes  the  statistical  interpretation  of  results, 
and,  for  this  reason,  comprehensive  printed  reports  of  auxiliary 
statistical  information,  often  in  graphical  form,  are  automatically 
provided  to  augment  the  basic  statistical  computations  performed  by 
each  user-callable  STARPAC  subroutine.  STARPAC  thus  provides  the 
best  features  of  many  stand-alone  statistical  software  programs 
within  the  flexible  environment  of  a  subroutine  library. 

This  Note  documents  16  subroutines  for  nonlinear  least  squares 
regression.  Twelve  of  these  compute  the  least  squares  estimates, 
performing  either  weighted  or  unweighted  analysis  with  either 
numerically  approximated  or  user-supplied  (analytic)  derivatives. 
The  other  four  are  user-callable  subroutines  for  two  procedures  used 
within  the  estimation  code:  the  first  selects  optimum  step  sizes 
for  approximating  the  partial  derviatives  of  the  model;  and  the 
second  checks  the  validity  of  a  user-supplied  derivative 
subroutine. 


Key  words:  derivative  checking;  NL2S0L;  nonlinear  least  squares; 
nonlinear  regression;  quasi-Newton  methods;  STARPAC;  statistical 
computing;  statistical  subroutine  library;  statistics;  derivative 
step  size  selection;  weighted  nonlinear  least  squares. 
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CHAPTER  1 
NONLINEAR  LEAST  SQUARES  REGRESSION 


A.   Introduction 

STARPAC  contains  16  user-callable  subroutines  for  nonlinear  least  squares 
regression.  Twelve  of  these  are  estimation  subroutines  that  compute  the  least 
squares  solution  as  described  below,  performing  either  weighted  or  unweighted 
regression  with  either  numerically  approximated  or  user-supplied  (analytic) 
derivatives.  The  estimation  subroutines  allow  three  levels  of  control  of  the 
computations  and  printed  output,  and  allow  the  user  to  specify  a  subset  of  the 
parameters  to  be  treated  as  constants,  with  their  values  held  fixed  at  their 
input  values.  This  last  feature  allows  the  user  to  examine  the  results 
obtained  estimating  various  subsets  of  the  parameters  of  a  general  model 
without  rewriting  the  model  subroutine  for  each  subset.  The  other  four 
subroutines  described  in  this  chapter  are  utility  procedures  which  choose 
optimum  step  sizes  for  numerically  approximating  the  derivative,  and  which 
verify  the  correctness  of  user-supplied  (analytic)  derivatives. 

Each  of  the  subroutines  described  in  this  chapter  assumes  that  the 
observations  of  the  dependent  variable,  y^ ,  are  modeled  by 

y±   =   fi(xi>£)  +  £±  for  1-1,  .-.,  N, 

where 

N       is  the  number  of  observations; 

f^      is   the  function  (nonlinear  in  its  parameters)  that  models   the  itn 
observation; 

x^      is  the  vector  of  the  M  independent  variables  at  the  itn  observation; 

3       is  the  vector  of  the  NPAR  model  parameters;  and 

e^  is   the  unobservable  random  error  in   the  itn  observation,   which  is 

estimated  by  the  itn  residual. 

The  least  squares  estimates  of  the  parameters,  6,  are  obtained  using  an 
iterative  procedure  that  requires  the  matrix  of  partial  derivatives  of  the 
model  with  respect  to  each  parameter, 

D(i,k)  =  3f1(x1,3)/33(k)     for  1=1,  ...,  N 

and  k  =  1,  . . . ,  NPAR. 

The     derivative       matrix     may       be     supplied        analytically        or     approximated 
numerically . 
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The  least  squares  solution,  3,  is  that  which  minimizes  (with  respect  to 
3)  the  residual  sum  of  squares  function, 

N  N 

RSS(3)   =     I       [{yi  ~   f±(xit&))2-wt±]    =     I       [e±2-wt±l 

i=l  1=1 

where  carat  (*)  denotes  the  estimated  quantity,  and 

wt^  is  the  weight  assigned  to  the  ifc"  observation  (wt^  =  1.0  in  the 
"unweighted"  case).  Appendix  A  discusses  several  common  applications 
for  weighted  least  squares. 

The  user  must  supply  both  initial  values  for  the  parameters,  and  a 
subroutine  [see  §D,  argument  NLSMDL]  to  compute  f^x^,^),  1  =  1,  .  .  . ,  N, 
i.e.,  the  predicted  values  of  the  dependent  variable  given  the  independent 
variables  and  the  parameter  values  from  iteration  I,  %  =  1,  2,  ...  .  Initial 
parameter  values  should  be  chosen  with  care,  since  good  values  can 
significantly  reduce  computer  time. 

STARPAC  provides  a  variety  of  subroutines  to  accommodate  many  levels  of 
user  sophistication  and  problem  difficulty.  Users  are  directed  to  §B  for  a 
brief  description  of  the  features  and  subroutines  available.  The  actual 
declaration  and  CALL  statements  are  given  in  §C,  and  the  subroutine  arguments 
are  defined  in  §D.  Generally,  the  computational  details  provided  in  §E  are 
not  needed  unless  difficulties  arise.  Sample  programs  and  their  output  are 
shown  in  %¥ . 


B .   Subroutine  Descriptions 

B .1   Nonlinear  Least  Squares  Estimation  Subroutines 

The  simplest  of  the  12  nonlinear  least  squares  estimation  subroutines, 
NLS,  requires  neither  user-supplied  weights  nor  analytic  derivatives.  The 
estimated  results  and  a  variety  of  statistics  are  automatically  summarized  in 
a  five-part  printed  report,  and  the  estimated  parameters  and  residuals  are 
returned  to  the  user  via  the  subroutine  argument  list  (level  one  control, 
described  below).  Most  nonlinear  least  squares  problems  can  be  solved  using 
NLS. 

The  other  11  estimation  subroutines  add  the  weighting,  derivative,  and 
level  two  and  three  control  features  both  singly  and  in  combination,  providing 
greater  flexibility  to  the  user  at  the  price  of  less  simplicity.  These 
features  are  indicated  by  the  suffix  letter(s)  on  the  subroutine  name  (e.g., 
NLSS  and  NLSWDC). 
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•  Suffix  W  indicates  user-supplied  weights. 

•  Suffix  D  indicates  user-supplied  (analytic)  derivatives. 

•  Suffix  C  indicates  level  two  control  of  the  computations. 

m   Suffix  S  indicates  level  three  control  of  the  computations. 

The  three  levels  of  computation  and  printed  output  control  are  described 
as  follows. 

o  In  level  one,  a  five-part  printed  report  is  automatically  provided 
[see  §E.l.b],  and  the  estimated  model  parameters  and  residuals  are 
returned  to  the  user  via  the  argument  list. 

•  Level  two  also  returns  the  estimated  parameters  and  residuals,  and, 
in  addition,  allows  the  user  to  supply  arguments  to  indicate 

-  a  subset  of  the  model  parameters  to  be  treated  as  constants, 
with  their  values  held  fixed  at  their  input  values; 

-  either  the  step  sizes  used  to  compute  the  numerical 
approximations  to  the  derivative  [see  §E.2],  or,  when 
user-supplied  analytic  derivatives  are  used,  whether  they  will 
be  checked  [see  §E.3]; 

-  the  maximum  number  of  iterations  allowed; 

-  the  convergence  criteria; 

-  the  scale  (i.e.,  the  typical  size)  of  each  parameter; 

-  the  maximum  change  allowed  in  the  parameters  at  the  first 
iteration; 

-  how  the  variance-covariance  matrix  is  to  be  approximated;  and 

-  the  amount  of  printed  output  desired. 

•  Level  three  has  all  the  features  of  level  two,  and,  in  addition 
returns  the  following  estimated  values  via  the  argument  list: 

-  the  number  of  nonzero  weighted  observations  (only  when  a 
weighted  analysis  is  performed); 

-  the  number  of  parameters  actually  estimated; 

-  the  residual  standard  deviation; 

-  the  predicted  values; 

-  the  standard  deviations  of  the  predicted  values; 

-  the  standardized  residuals;  and 

-  the  variance-covariance  matrix  of  the  estimated  parameters. 

B .2   Derivative  Step  Size  Selection  Subroutines 

When  the  partial  derivatives  used  in  the  nonlinear  least  squares  solution 
are  not  available  analytically,  STARPAC  subroutines  approximate  them 
numerically.  In  this  case,  the  subroutines  can  select  optimum  step  sizes  for 
approximating  the  derivatives  [see  §E.2].  The  user  also  has  the  option  of 
computing  these  step  sizes  independently  of  the  estimation  process  by  calling 
either  of  the  two  step  size  selection  subroutines  directly.   For  example,  when 
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planning  to  use  the  parameter  fixing  capability  [see  §D,  argument  IFIXED]  to 
examine  several  subsets  of  the  parameters  of  a  general  model,  computing  the 
step  sizes  first  and  passing  them  to  the  estimation  subroutine  is  more 
efficient  than  recomputing  them  each  time  the  estimation  subroutine  is 
called. 

The  simplest  of  the  two  user-callable  step  size  selection  subroutines, 
STPLS,  summarizes  the  step  size  selection  information  for  each  parameter  in  a 
printed  report,  and  returns  the  step  sizes  to  the  user  via  the  subroutine 
argument  list. 

The  second  step  size  selection  subroutine,  STPLSC,  differs  from  STPLS 
only  in  that  it  enables  the  user  to  supply  arguments  to  indicate 

-  the  number  of  reliable  digits  in  the  model  results; 

-  the  number  of  exemptions  allowed  by  the  acceptance  criteria; 

-  the  scale  (i.e.,  the  typical  size)  of  each  parameter;  and 

-  the  amount  of  printed  output  desired. 


B .3   Derivative  Checking  Subroutines 

When  the  partial  derivatives  used  in  the  nonlinear  least  squares  solution 
are  available  analytically,  the  user  can  code  them  for  use  by  the  estimation 
subroutines  [see  §D,  argument  NLSDRV] .  Because  coding  errors  are  a  common 
problem  with  user-supplied  derivatives,  the  STARPAC  estimation  subroutines 
automatically  check  the  validity  of  the  user-supplied  derivative  code  by 
comparing  its  results  to  numerically  approximated  values  for  the  derivative. 
When  the  results  are  questionable,  the  checking  procedure  attempts  to 
determine  whether  the  problem  lies  with  the  user's  code  or  with  the  accuracy 
of  the  numerical  approximation  [see  §E.3].  Although  the  checking  procedure  is 
automatically  available  to  the  estimation  subroutines  which  accept 
user-supplied  derivatives,  the  user  may  want  to  check  the  derivative  code 
independently  of  the  estimation  process.  In  these  cases,  the  user  can  call 
either  of  the  two  derivative  checking  subroutines  directly,  and  suppress 
checking  by  the  estimation  subroutines.   [See  §D,  argument  IDRVCK. ] 

The   simplest   of   the  two   derivative   checking   subroutines,   DCKLS, 
summarizes  the  results  of  the  check  in  a  printed  report. 

The  second  derivative  checking  subroutine,  DCKLSC,   differs   from  DCKLS 

only  in  that  it  enables  the  user  to  supply  arguments  to  indicate 

-  the  number  of  reliable  digits  in  the  model  results; 

-  the  agreement  tolerance; 

-  the  scale  (i.e.  ,  the  typical  size)  of  each  parameter; 

-  the  row  at  which  the  derivative  is  to  be  checked;  and 

-  the  amount  of  printed  output  desired. 
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C .   Subroutine  Declaration  and  CALL  Statements 

NOTE:  Argument  definitions  and  sample  programs  are  given  in  §D  and  §F , 
respectively.  The  conventions  used  to  present  the  following  declaration  and 
CALL  statments  are  given  in  the  Introduction  to  STARPAC  [Donaldson  and  Tryon, 
1983,  chapter  1,  §B,  and  §D]  . 


Nonlinear  Least  Squares  Estimation  Subroutines 

The  <basic  declaration  block>  identifies  declaration  statements  that  are 
needed  by  all  of  the  nonlinear  least  squares  estimation  subroutines.  The  user 
should  substitute  the  following  four  statements  for  each  occurrence  of  <basic 
declaration  block>  given  below. 

<real>  Y(n),  XM(n,r?),  ?AR(npar) ,  RES(^?) 
DOUBLE  PRECISION  DSTAK(  Idstak) 
COMMON  /CSTAK/  DSTAK 
EXTERNAL  NLSMDL 


NLS:         Unweighted     nonlinear  least  squares  regression;  numerically 

approximated  derivatives;  five-part  printed  report;  minimal  return  of 
results. 

<basic  declaration  block> 


CALL  NLS  (Y,  XM ,  N,  M,  IXM ,  NLSMDL,  PAR,  NPAR,  RES,  LDSTAK) 


NLSC:       Unweighted     nonlinear  least  squares  regression;  numerically 

approximated  derivatives;  user-supplied  control  values;  minimal  return 
of  results. 

<basic  declaration  block> 

INTEGER  IFIXED( npar) 

<real>  ST? (n par) ,    STOPSS,  STOPP,  SCALE(npar) ,  DELTA 


CALL  NLSC  (Y,  XM ,  N,  M,  IXM,  NLSMDL,  PAR,  NPAR,  RES,  LDSTAK, 
1    IFIXED,  STP,  MIT,  STOPSS,  STOPP,  SCALE,  DELTA,  IVAPRX,  NPRT) 
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NLSS:       Unweighted     nonlinear  least  squares  regression;  numerically 

approximated  derivatives;  user-supplied  control  values;  maximal  return 
of  results. 

<basic  declaration  block> 

INTEGER  IFIXEV(  npar) 

<real>    ST? (n par) ,    STOPSS,    STOPP,    SCKLE(npar) ,    DELTA 

<real>  RSD,    PV(  n)  ,   SDPV(  n)  ,   SDRES(w),   VCV(npare3npare) 


CALL  NLSS    (Y,    XM ,    N,    M,    IXM ,    NLSMDL,    PAR,    NPAR,    RES,    LDSTAK, 

1  IFIXED,    STP,    MIT,    STOPSS,    STOPP,    SCALE,    DELTA,    IVAPRX,    NPRT, 

2  NPARE,    RSD,    PV,    SDPV,    SDRES,    VCV,    IVCV) 


NLSW:       Weighted     nonlinear  least  squares  regression;  numerically  approximated 
derivatives;  five-part  printed  report;  minimal  return  of  results. 

<basic  declaration  block> 
<real>  WT(rz) 


CALL  NLSW  (Y,  WT,  XM ,  N,  M,  IXM,  NLSMDL,  PAR,  NPAR,  RES,  LDSTAK) 


NLSWC:     Weighted     nonlinear  least  squares  regression;  numerically  approximated 
derivatives;  user- supplied,  control  values;  minimal  return  of  results. 

<basic  declaration  block> 

INTEGER  IFIXED(rzpar) 

<real>  WT(rz) 

<real>  STP (n par)  ,  STOPSS,  STOPP,  SCALE( npar)  ,   DELTA 


CALL  NLSWC  (Y,  WT,  XM ,  N,  M,  IXM,  NLSMDL,  PAR,  NPAR,  RES,  LDSTAK, 
1    IFIXED,  STP,  MIT,  STOPSS,  STOPP,  SCALE,  DELTA,  IVAPRX,  NPRT) 
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NLSWS:     Weighted     nonlinear  least  squares  regression;  numerically  approximated 
derivatives;  user-supplied  control  values;  maximal  return  of  results, 

<basic  declaration  block> 

INTEGER  IFIXED(  n par) 

<real>  WT(n) 

<real>    STP(npar),    STOPSS,    STOPP,    SCALE(npar) ,   DELTA 

<real>  RSD,    PV(  n)  ,   SDPV(  n) ,   SDRES(n),   VCV  ( n par esn pare) 


CALL  NLSWS    (Y,    WT,    XM ,   N,    M,    IXM ,    NLSMDL,    PAR,    NPAR,   RES,    LDSTAK, 

1  IFIXED,    STP,    MIT,    STOPSS,    STOPP,    SCALE,    DELTA,    IVAPRX,    NPRT, 

2  NNZW,   NPARE,    RSD,    PV,    SDPV,    SDRES,    VCV,    I VCV) 


NLSD:        Unweighted     nonlinear  least  squares  regression;  user-supplied 

derivatives;  five-part  printed  report;  minimal  return  of  results, 

<basic  declaration  block> 
EXTERNAL  NLSDRV 


CALL  NLSD  (Y,  XM,  N,  M,  IXM,  NLSMDL,  NLSDRV,  PAR,  NPAR,  RES, 
1    LDSTAK) 


NLSDC:     Unweighted     nonlinear  least  squares  regression;  user-supplied 

derivatives;  user-supplied  control  values;  minimal  return  of  results, 

<basic  declaration  block> 

EXTERNAL  NLSDRV 

INTEGER  IFIXED(  npar) 

<real>  STOPSS,  STOPP,  SCALE(npar),  DELTA 


CALL  NLSDC  (Y,  XM ,  N,  M,  IXM,  NLSMDL,  NLSDRV,  PAR,  NPAR,  RES, 

1  LDSTAK,  IFIXED,  IDRVCK,  MIT,  STOPSS,  STOPP,  SCALE,  DELTA, 

2  IVAPRX,  NPRT) 
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NLSDS:     Unweighted     nonlinear  least  squares  regression;  user-supplied 

derivatives;  user-supplied  control  values;  maximal  return  of  results, 

<basic  declaration  block> 

EXTERNAL  NLSDRV 

INTEGER  IFIXED(npa^) 

<real>  STOPSS,  STOPP,  SCALE(npar),  DELTA 

<real>  RSD,  PV(w),  SDPV(  n) ,  SDRES(n),  VCV ( n par e3n pare) 


CALL  NLSDS  (Y,  XM ,  N,  M,  IXM ,  NLSMDL,  NLSDRV,  PAR,  NPAR,  RES, 

1  LDSTAK,  IFIXED,  TDRVCK,  MIT,  STOPSS,  STOPP,  SCALE,  DELTA, 

2  IVAPRX,  NPRT,  NPARE,  RSD,  PV,  SDPV,  SDRES,  VCV,  I VCV) 


NLSUD:      Weighted     nonlinear  least  squares  regression;  user-supplied 

derivatives;  five-part  printed  report;  minimal  return  of  results. 

<basic  declaration  block> 
EXTERNAL  NLSDRV 
<real>  WT(rz) 


CALL  NLSWD  (Y,  WT,  XM ,  N,  M,  IXM,  NLSMDL,  NLSDRV,  PAR,  NPAR,  RES, 
1    LDSTAK) 


NLSWDC:  Weighted     nonlinear  least  squares  regression;  user-supplied 

derivatives;  user-supplied  control  values;  minimal  return  of  results, 

<basic  declaration  block> 

EXTERNAL  NLSDRV 

INTEGER  IFIXED(  npar) 

<real>  WTO) 

<real>  STOPSS,  STOPP,  SCALE(npar),    DELTA 


CALL  NLSWDC  (Y,  WT,  XM ,  N,  M,  IXM,  NLSMDL,  NLSDRV,  PAR,  NPAR,  RES, 

1  LDSTAK,  IFIXED,  IDRVCK,  MIT,  STOPSS,  STOPP,  SCALE,  DELTA, 

2  IVAPRX,  NPRT) 
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NLSWDS:  Weighted     nonlinear  least  squares  regression;  user-supplied 

derivatives;  user-supplied  control  values;  maximal  return  of  results, 

<basic  declaration  block> 

EXTERNAL  NLSDRV 

INTEGER  IFIXED(npar) 

<real>   WT(w) 

<real>  STOPSS,  STOPP,  SCALE(  npar) ,  DELTA 

<real>  RSD,  PV(  n) ,  SDPV(  n) ,  SDRES(n),  VCV (n pare }n pare) 


CALL  NLSWDS  (Y,  WT,  XM ,  N,  M,  IXM,  NLSMDL,  NLSDRV,  PAR,  NPAR,  RES, 

1  LDSTAK,  IFIXED,  IDRVCK,  MIT,  STOPSS,  STOPP,  SCALE,  DELTA, 

2  IVAPRX,  NPRT,  NNZW,  NPARE,  RSD,  PV,  SDPV,  SDRES,  VCV,  IVCV) 


Step  Size  Selection  Subroutines 

STPLS:     Optimum     step  size  selection  for  numerically  approximating 
derivatives;  printed  report;  maximal  return  of  results. 

<real>  XM(n,w),  FkR(npar) ,    ST?(npar) 
DOUBLE  PRECISION  DSTAK(  Idstak) 
COMMON  /CSTAK/  DSTAK 
EXTERNAL  NLSMDL 


CALL  STPLS  (XM,  N,  M,  IXM,  NLSMDL,  PAR,  NPAR,  LDSTAK,  STP) 


S^PLSC:  Optimum     step  size  selection  for  numerically  approximating 

derivatives;  user-supplied  control  values;  maximal  return  of  results, 

<real>  m(n,m)  ,    PAR( npar)  y    STP(npar) 
<real>  EXMPT,  SCALE(npar) 
DOUBLE  PRECISION  DSTAK(  Idstak) 
COMMON  /CSTAK/  DSTAK 
EXTERNAL  NLSMDL 


CALL  STPLSC  (XM,  N,  M,  IXM,  NLSMDL,  PAR,  NPAR,  LDSTAK,  STP, 
1    NETA,  EXMPT,  SCALE,  NPRT) 
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Derivative  Checking  Subroutines 

DCKLS:      Derivative  checking;  -printed  report;  maximal  return  of  results, 

<real>  XM(n,m),  PAR( npar)  t    ST?(npar) 
DOUBLE  PRECISION  DSTAK(  Idstak) 
COMMON  /CSTAK/  DSTAK 
EXTERNAL  NLSMDL,  NLSDRV 

CALL  DCKLS  (XM,  N,  M,  IXM ,  NLSMDL,  NLSDRV,  PAR,  NPAR,  LDSTAK) 


DCKLSC:  Derivative  checking;  user-supplied  control  values;  maximal  return  of 
results. 

<real>  XM(n,m),  PAR( npar),    STP(npar) 

<real>  SCALE(npar) 

DOUBLE  PRECISION  DSTAK(  Idstak) 

COMMON  /CSTAK/  DSTAK 

EXTERNAL  NLSMDL,  NLSDRV 


CALL  DCKLSC  (XM ,  N,  M,  IXM,  NLSMDL,  NLSDRV,  PAR,  NPAR,  LDSTAK, 
1    NETA,  NTAU,  SCALE,  NROW ,  NPRT) 


D.   Dictionary  of  Subroutine  Arguments  and  COMMON  Variables 

NOTE:    — >  indicates   that  the  argument  is  input  to  the  subroutine,  and  that 

the  input  value  is  preserved; 
< —  indicates  that  the  argument  is  returned  by  the  subroutine; 
<->  indicates  that  the  argument   is  input   to  the  subroutine   and  that 

the  input  value  will  be  overwritten  by  the  subroutine; 
indicates  that  the  argument   is  input   to  some  subroutines   and  is 

returned  by  others; 
***  indicates  that  the  argument  is  a  subroutine  name; 
•••  indicates  that  the  variable  is  passed  via  COMMON. 


< —  The  matrix  of  exact  dimension  N  by  NPAR  that  contains  the  partial 
derivatives  of  the  model  with  respect  to  each  parameter, 
PAR(k) ,  k  =  1,  ...,  NPAR.  This  argument  is  used  within  derivative 
subroutine  NLSDRV  [see  argument  NLSDRV  below]. 
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DELTA  — >  The  maximum  scaled  change  allowed  In  the  parameters  at  the  first 
iteration,  i.e.,  6q  [see  §E.l.a].  The  default  value  is  100.0. 
When  DELTA  <  0.0,  or  when  DELTA  is  not  an  argument  of  the 
subroutine  CALL  statement,  the  default  value  is  used.  A  smaller 
value  of  DELTA  may  be  appropriate  if,  at  the  first  iteration,  the 
computation  of  the  predicted  values  from  the  user's  model 
subroutine  produces  an  arithmetric  overflow,  or  the  parameters 
leave  the  region  of  interest  in  parameter  space.  A  reasonable 
alternative  to  the  default  value  of  DELTA  is  an  upper  bound  to  the 
scaled  change  that  the  estimated  parameters  should  be  allowed  to 
make  on  the  first  iteration, 

DELTA  =  min{ |AmaxB(k) |/SCALE(k),  for  k  =  1,  ...,  NPAR} 

where  Amax3(k)  is  the  maximum  change  allowed  for  the  ktn  parameter 
at  the  first  iteration. 


DSTAK 


The  DOUBLE  PRECISION  vector  in  COMMON  /CSTAK/  of  dimension  at 
least  LDSTAK.  DSTAK  provides  workspace  for  the  computations.  The 
first  LDSTAK  locations  of  DSTAK  will  be  overwritten  during 
subroutine  execution. 


EXMPT  — >  The  proportion  used  to  compute  the  number  of  observations, 
a  =  EXMPT»N,  for  which  the  forward  difference  quotient  derivative 
with  respect  to  a  given  parameter  is  exempted  from  meeting  the 
acceptance  criteria  for  step  size  selection  [see  §E.2.a].  The 
default  value  for  EXMPT  is  0.1  (10  percent).  When  the 
user-supplied  value  is  outside  the  range  [0.0,  1.0],  or  when  EXMPT 
is  not  an  argument  of  the  subroutine  CALL  statement,  the  default 
value  is  used. 

IDRVCK  — >  The  indicator  variable  used  to  designate  whether  or  not  the 
user-supplied  derivative  subroutine  is  to  be  checked.  When 
IDRVCK  *  0  the  derivative  is  checked,  and  when  IDRVCK  =0  it  is 
not.  The  default  value  is  IDRVCK  *  0.  When  IDRVCK  is  not  an 
argument  of  the  subroutine  CALL  statement,  the  default  value  is 
used. 


IERR 


An  error  flag  returned  in  COMMON  /ERRCHK/.  [See  Donaldson  and 
Tryon,  1983,  chapter  1,  §D.5.]  Note  that  using  (or  not  using)  the 
error  flag  will  not  affect  any  error  messages  that  are 
automatically  printed,  even  when  the  user  has  suppressed  the 
normal  printed  output. 


For  the  estimation  subroutines: 

IERR  =  0  indicates  that  no  errors  were  detected,  and   that  the 

iterations  converged  satisfactorily. 
IERR  =  1  indicates   that   improper   input  was   detected. 

continued  — 
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IERR  =  2  indicates  that  the  computation  of  the  residual  sum  of 
squares  using  the  initial  parameter  values  produced  an 
arithmetic  overflow.  The  user  should  reduce  the  size 
of  DELTA  or  should  supply  new  starting  values. 

IERR  =  3  indicates  that  the  model  is  computationally  singular, 
which  means  the  model  has  too  many  parameters  near  the 
solution.  The  user  should  examine  the  model  and  data 
to  determine  and  remove  the  cause  of  the  singularity. 

IERR  =  4  indicates  that  at  least  one  of  the  standardized 
residuals  could  not  be  computed  because  its  standard 
deviation  was  zero.  The  validity  of  the  covariance 
matrix  is  questionable. 

IERR  =  5  indicates  false  convergence.   [See  §E.l.a.] 

IERR  =  6  indicates  that  convergence  was  not  reached  in  the 
allowed  number  of  iterations  or  model  subroutine 
calls.   [See  argument  MIT.] 

IERR  =  7  indicates  that  the  variance-covariance  matrix  could 
not  be  computed. 

For  the  step  size  selection  subroutines: 

IERR  =  0  indicates  that  no  errors  were  detected,  and  that  all 
the  step  sizes  satisfied  the  selection  criteria. 

IERR  =  1  indicates   that   improper   input  was   detected. 

IERR  =  2  indicates  that  one  or  more  of  the  step  sizes  did  not 
satisfy  the  selection  criteria. 

For  the  derivative  checking  subroutines: 

IERR  =  0  indicates  that  no  errors  were  detected,  and  that  the 
user-supplied  derivative  code  appears  to  be  correct. 

IERR  =  1  indicates   that  improper   input  was   detected. 

IERR  =  2  indicates  that  the  user-supplied  derivative  code  and 
numerical  derivatives  do  not  agree  for  at  least  one 
parameter,  but  that  in  each  case  of  disagreement  the 
accuracy  of  the  numerical  derivatives  is  questionable. 
Further  testing  is  suggested. 

IERR  =  3  indicates  that  the  user-supplied  derivative  code  and 
numerical  derivatives  do  not  agree  for  at  least  one 
parameter,  and  in  at  least  one  instance  of 
disagreement  there  is  no  reason  to  doubt  the  numerical 
derivatives. 

IFIXED  — >  The  vector  of  dimension  at  least  NPAR  that  contains  values  used  to 
indicate  whether  the  corresponding  parameter  in  PAR  is  to  be 
treated  as  a  fixed  constant  or  is  to  be  estimated.  If 
IFIXED(I)  >  0,  PAR(I)  will  be  held  fixed  at  its  input  value;  if 
IFIXED(I)  =  0,  PAR(I)  will  be   estimated   using  the   least  squares 

continued  — 
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IVCV 


procedure  described  in  §A.  The  default  values  are  IFIXED(I)  =  0, 

I  =   1,  •••,  NPAR,   i.e.,  all   parameters   are   estimated.    When 

IFIXED(l)  <  -1,   or  when  IFIXED   is   not   an  argument   of   the 

subroutine  CALL  statement,  the  default  value  will  be  used. 

•>  The  exact   value  of   the   first  dimension  of  the  matrix  VCV   as 
specified  in  the  calling  program. 


IVAPRX  — >  The  indicator  variable  used  to  specify  how  the  variance-covariance 
matrix,  VCV,  is  to  be  approximated.  Three  approximations  are 
available: 

(1)  VCV    =   RSD2«(DT.W.D) 

(2)  VCV    =   RSD2.H~1 

(3)  VCV    =   RSD2.H_1.(DT.W.D).H-1 
where 

H  =  DT.W.D  +  {  y   ei'wt.  .(82ei/3B(j)96(k))  for  j  =  1,  ...,  NPAR 
i=l  and  k  =  1,  . . . ,  NPAR} ; 

W   is   an  N  by   N  diagonal  matrix  of   user-supplied  weights, 

W  =  diag{wtis  i  =  1,  ...,  N}, 

when  a  weighted  analysis  is  performed,  and  is  the  identity  matrix 
otherwise. 


Approximation  (1)  is  based  on  the  assumption  that  H  ~  D  »W«D 
because  the  residuals  are  sufficiently  small  at  the  solution; 
approximation  (2)  is  based  on  the  assumption  that  the  necessary 
conditions  for  asymptotic  maximum  likelihood  theory  have  been  met; 
and  approximation  (3)  is  based  on  the  assumption  that  the 
necessary  conditions  for  asymptotic  maximum  likelihood  theory  may 
be  violated.  All  approximations  to  the  variance-covariance  matrix 
are  subject  to  sampling  variation,  however,  because  they  are 
computed  using  the  estimated  parameter  values.  The 
variance-covariance  matrix  computed  for  any  particular  nonlinear 
least  squares  solution  should  thus  be  regarded  as  only  a  rough 
estimate,  correct  to  within  an  order  of  magnitude  [Bard,  1974]. 

If  IVAPRX  =  1  or  4  then  approximation  (1)  is  used; 

=  2  or  5  then  approximation  (2)  is  used;  and 
=  3  or  6  then  approximation  (3)  is  used. 

If  IVAPRX  =1,  2,  or  3,  then,  when  user-supplied  analytic 
derivatives  are  available  [see  argument  NLSDRV] ,  they  are  used  to 
compute  VCV;  if  IVAPRX  =4,  5,  or  6,  then  only  the  predicted 
values  from  the  model  subroutine  are  used  to  compute  VCV.  When 
analytic  derivatives  are  available,  options  1,  2,  or  3,  will 
generally  result  in  a  faster,  more  accurate  computation  of  VCV. 

continued  — 
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The  default  value  for  IVAPRX  is  1.  When  argument  IVAPRX  is 
outside  the  range  [1,  6],  or  when  IVAPRX  is  not  an  argument  of  the 
subroutine  CALL  statement,  then  the  default  value  will  be  used. 

IXM  — >  The  exact  value  of  the  first  dimension  of  the  matrix  XM  as 
specified  in  the  calling  program. 

LDSTAK  — >  The  length  of  the  DOUBLE  PRECISION  workspace  vector  DSTAK.  LDSTAK 
must  equal  or  exceed  the  appropriate  value  given  below.  If  the 
single  precision  version  of  STARPAC  is  being  used,  P  =  0.5 
otherwise  P  =  1.0  [see  Donaldson  and  Tryon,  1983,  chapter  1,  $B] . 

For  NLS,  NLSC,  NLSS,  NLSW,  NLSWC,  NLSWS: 

LDSTAK  >   27  +  max) IS* (N+NPAR) ,  30+NPAREJ  + 

P.max{lS'10.N,  94+N« (3+NPAR)+(3 .NPARE2+37.NPARE) /2 } 

where  IS  =  1  if  default  values  are  used  for  the  derivative 
step  sizes,  and  IS  =  0  otherwise. 

For  NLSD,  NLSDC ,  NLSDS,  NLSWD,  NLSWDC ,  NLSWDS: 

LDSTAK  >  45  +  NPAR  +  P«  (94+N  •  (3+NPAR)+(3  .NPARE2+35»  NPARE)  /2  ) 

For  STPLS,  STPLSC: 

LDSTAK  >  27  +  (N+NPAR)  +  P«10-N 

For  DCKLS,  DCKLSC: 

LDSTAK  >    14  +  NPAR  +  P-  (N-NPAR+N+NPAR) 

M  — >  The  number  of  independent  variables,  i.e.,  the  number  of  columns 
of  data  in  XM. 

MIT  — >  The  maximum  number  of  iterations  allowed.  This  argument  is  also 
used  to  compute  the  maximum  number  of  model  subroutine  calls, 
(2»MIT).  The  iterations  will  stop  if  either  limit  is  reached, 
although,  as  a  rule,  the  maximum  number  of  iterations  will  be 
reached  first.  The  default  value  for  the  maximum  number  of 
iterations  is  21.  When  MIT  <  0,  or  when  MIT  is  not  an  argument  of 
the  subroutine  CALL  statement,  the  default  value  will  be  used. 

N       — >  The  number  of  observations. 
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NETA  — >  The  number  of  reliable  decimal  digits  in  the  predicted  values  (PV) 
computed  by  the  user's  model  subroutine.  The  default  value  for 
NETA  is  experimentally  determined  by  the  procedure  described  in 
Appendix  B.  The  default  value  will  be  used  when  NETA  is  not  an 
argument  in  the  subroutine  CALL  statement,  or  when  the 
user-supplied  value  of  NETA  is  outside  the  range  [1,  DIGITS], 
where  DIGITS  is  the  number  of  decimal  digits  carried  by  the  user's 
computer  for  a  single  precision  value  when  the  single  precision 
version  of  STARPAC  is  being  used  and  is  the  number  carried  for  a 
double  precision  value  otherwise. 

NLSDRV  ***  The  name  of  the  user-supplied  subroutine  that  computes  the  partial 
derivative  matrix  (Jacobian) .  This  argument  must  be  listed  in  an 
EXTERNAL  statement  in  the  program  which  calls  the  STARPAC 
estimation  or  derivative  checking  subroutine.  The  form  of  the 
derivative  subroutine  argument  list  and  dimensioning  statements 
must  be  exactly  as  shown  below,  although  if  there  is  only  one 
independent  variable  (M  =  1),  XM  may  be  declared  to  be  a  vector 
with  dimension  IXM . 

SUBROUTINE  NLSDRV  (PAR,  NPAR,  XM ,  N,  M,  IXM,  D) 
<real>  PAR(NPAR),  XM(IXM,M),  D(N,NPAR) 

<  Computations  for  D(I,J), 

1=1,  . . . ,  N  and  J  =  1 ,  . . . ,  NPAR  > 

RETURN 
END 

NLSMDL  ***  The  name  of  the  user-supplied  subroutine  that  computes  the 
predicted  value  of  the  dependent  variable  given  the  independent 
variables  and  the  current  values  of  the  model  parameters.  This 
argument  must  be  listed  in  an  EXTERNAL  statement  in  the  program 
which  calls  the  STARPAC  estimation,  step  size  selection,  and/or 
derivative  checking  subroutines.  The  form  of  the  model  subroutine 
argument  list  and  dimensioning  statements  must  be  exactly  as  shown 
below,  although  if  there  is  only  one  independent  variable  (M  =  1) , 
XM  may  be  declared  to  be  a  vector  with  dimension  IXM. 

SUBROUTINE  NLSMDL  (PAR,  NPAR,  XM ,  N,  M,  IXM,  PV) 
<real>  PAR(NPAR) ,  XM(IXM,M),  PV(N) 

<  Computations  for  PV(I),  1=1,  ...,  N  > 

RETURN 
END 

NNZW    < —  The  number  of  observations  with  nonzero  weights. 
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NPAR  — >  The  number  of  parameters  in  the  model,  including  both  those  held 
fixed  at  their  starting  values,  and  those  which  are  to  be 
estimated. 

NPARE  < —  The  number  of  parameters  actually  estimated,  i.e.,  the  number  of 
zero  elements  in  IFIXED.  N.B.  This  value  is  returned  by  the 
estimation  subroutines. 

NPRT    — >  The  argument  controlling  printed  output. 

For  the  estimation  subroutines: 

NPRT  is   a  five-digit  integer,   in  which   the  value  of   the  1^ 
digit  (counting  from  left  to  right)  is  used  to  control   the  Itn 
section  of  the  output. 

If  the  Itn  digit  =  0  the  output  from  the  Itn  section  is 

suppressed; 
=  1  the  brief  form  of  the  Itn  section  is 

given; 
>   2  the  full  form  of  the  Itn  section  is 
given. 
The  default   value  for  NPRT  is  11112.    When  NPRT  <  -1,  or  when 
NPRT  is  not  an  argument  in  the  subroutine  CALL  statement,  the 
default  value  will  be  used.    If  the  convergence  criteria  are 
not  satisfied,   the  subroutine  gives  a   suitable  warning  and 
provides  a  printed  report  even  if  NPRT  =0.    A  full  discussion 
of  the  printed  output  is  given  in  §E.l.b,  and  is   summarized  as 
follows. 

Section  1  lists  the   starting  estimates   and   control  values. 
Brief  output  and   full  output  are  the  same   for  this 
section. 
Section  2  reports  the  results  of  the  iterations.    Brief  output 
includes   information  only   about  the  first   and  last 
iteration  while   full   output   includes   information 
about  all  of  the  iterations. 
Section  3  provides  information  for  each  observation  based  on 
the   final   solution.     Brief   output   includes 
information  for  the  first  40  observations,  while  full 
output  provides  the  information  for  all  of  the  data. 
Section  4  is  a   set  of  four  residual  plots.    Brief  output  and 

full  output  are  the  same  for  this  section. 
Section  5  is   the  final   summary  of   the  estimated  parameters. 
Brief   output    does  not    include   printing   the 
variance-covariance  matrix,  while  full  output  does. 

For  the  step  size  selection  and  derivative  checking  subroutines: 

If  NPRT  =  0  the  output  is  suppressed. 
If  NPRT  *   0  the  output  is  provided. 

If  the  acceptance  criteria  were  not  met,  then  a   printed  report 
is  provided  even  if  NPRT  =  0. 
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NROW  — >  The  row  of  the  independent  variable  matrix  at  which  the 
user-supplied  derivative  code  is  to  be  checked.  The  default  value 
is  the  first  row  with  no  independent  variables  equal  to  zero;  when 
all  rows  have  one  or  more  independent  variables  equal  to  zero,  row 
one  will  be  used  for  the  default  value.  When  the  user-supplied 
value  is  outside  the  range  [1,  N]  ,  or  when  NROW  is  not  an  argument 
of  the  subroutine  CALL  statement,  the  default  value  will  be  used. 

NTAU  — >  The  agreement  tolerance,  i.e.,  the  number  of  digits  of  agreement 
required  between  the  user-supplied  derivatives  and  the  derivatives 
numerically  approximated  by  the  derivative  checking  subroutine. 
The  default  value  is  NETA/4 .  When  the  user-supplied  value  of  NTAU 
is  outside  the  range  [1,  NETA/2],  or  when  NTAU  is  not  an  argument 
of  the  subroutine  CALL  statement,  the  default  value  will  be"  tised. 


PAR 


The  vector  of  dimension  at  least  NPAR  that  contains  the  parameter 
values.  For  all  estimation  subroutines  it  must  contain  initial 
values  for  the  parameters  on  input  and  will  contain  the  final 
values  on  return.  For  the  step  size  and  derivative  checking 
subroutines  it  must  contain  the  parameter  values  at  which  the 
operations  are  to  be  performed. 


PV      < —  The   vector  of  dimension  at  least  N  that  contains   the  predicted 
values  of  the  dependent  variable  at  the  solution, 


PV(i)  =  fiCxi.3)  =  y±, 


for  i  =  1 ,  . . ,  N. 


RES  < —  The  vector  of  dimension  at  least  N  that  contains  the  residuals  at 
the  solution, 

RES(i)  =  y(i)  -  fiCxj.B)  -  y(i)  -  y(i)  =  e(i)  ,  for  i  =  1,  ...,  N. 

RSD     < —  The  residual  standard  deviation  at  the  solution, 

RSD  =  (RSS(3)/(NNZW-NPARE))1/2. 

SCALE  — >  The  vector  of  dimension  at  least  NPAR  that  contains  the  scale,  or 
typical  size,  of  each  parameter.  The  vector  SCALE  is  used  to 
normalize  the  size  of  each  parameter  so  that 


3(j)/SCALE(j) 


B(k)/SCALE(k) 


for  k  =  1,  . . .,  NPAR 
and  .1  =  1,  •  •  .  ,  NPAR. 


Values  of  |SCALE(k)|  >  I B(k) I  can  be  used  to  increase  the  step 
size  in  cases  where  the  model  function  is  known  to  be  insensitive 
to  small  changes  in  the  value  (3(k)  . 

For  the  estimation  subroutines: 

The  default  values  for  SCALE  are  selected  by  the  NL2S0L 
algorithm  [Dennis  et  al. ,  1981a  and  1981b]  and  are  updated  at 
each  iteration.     When  SCALE   is   not  an  argument   in  the 
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subroutine  CALL  statement,  or  when  the  user-supplied  value  for 
SCALE(l)  <  0,  the  default  procedure  will  be  used  to  select 
scale  values.  When  SCALE(l)  >  0,  values  of  SCALE(k)  <  0, 
k  =  2,  ...,  NPAR,  will  be  interpreted  as  an  input  error. 
User-supplied  scale  values  may  be  either  a  vector  of  the 
typical  size  of  each  parameter  or  a  vector  of  ones  if  the 
typical  sizes  of  the  parameters  are  roughly  equal; 
user-supplied  scale  values  can  sometimes  result  in  reduced 
computing  time  since  these  values  are  not  updated  at  each 
iteration. 

For  the  derivative  checking  and  step  size  selection  subroutines: 

The  default  values  of  SCALE  are  defined  as: 


SCALE(k)  =  1.0 


if  3(k)  =  0.0 


SCALE(k)  =  l3(k)|  otherwise 


for  k  =  1,  . .. ,  NPAR. 


When  SCALE  is  not  an  argument  in  the  subroutine  CALL  statement, 
or  when  the  user-supplied  value  of  |SCALE(k)|  <  l3(k)|,  the 
default  value  for  SCALE(k)  is  used.  When  SCALE(l)  <  0,  the 
default  values   will  be   used   for   SCALE(k),  k  =  1,  ...,  NPAR. 


When  SCALE(l)  >  0,  values   of   SCALE(k) 
will  be  interpreted  as  an  input  error. 


0,  k  =  2, 


NPAR, 


SDPV 


SDRES 


<— 


The  vector  of  dimension  at  least  N  that  contains  an  approximation 
to  the  standard  deviation  of  each  predicted  value  at  the 
solution, 


SDPV(i)  =  the  ith  diagonal  element  of  [rSD2. (D. VCV.Dr) ] 


1/2 


for  i  =  1,  ...,  N.  This  approximation  is  based  on  a  linearization 
of  the  model  in  the  neighborhood  of  the  solution;  the  validity  of 
the  approximation  depends  on  the  nonlinearity  of  the  model.  This 
approximation  may  be  extremely  inaccurate  for  a  problem  with  a 
highly  nonlinear  model. 

<--  The  vector  of  dimension  at  least  N  that  contains  an  approximation 
to  the  standardized  residuals  at  the  solution, 


SDRES(i)  =  RES(i)/[(RSD2/WT(i))  -  SDPV(i)2] 


1/2 


for  i  =  1,  .  .  .  ,  N,  which  is  the  itn  residual  divided  by  its 
individual  estimated  standard  deviation.  This  approximation  is 
based  on  a  linearization  of  the  model  in  the  neighborhood  of  the 
solution;  the  validity  of  the  approximation  depends  on  the 
nonlinearity  of  the  model.  This  approximation  may  be  extremely 
inaccurate  for  a  problem  with  a  highly  nonlinear  model. 
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STOPP  — >  The  stopping  value  for  the  convergence  test  based  on  the  maximum 
scaled  relative  change  in  the  parameters  at  the  most  recent 
iteration.  The  convergence  criterion  is  satisfied  if  the  current 
step  is  a  Newton  step  [see  Dennis  et  al. ,  1981a]  and 

max{|g8.1.1(k)-3g(k)|/SCALE(k)  fork  =  1,  ...,  NPAR}   <  ST0PP 
max{(le£+1(k)|+|B£(k)|  J/SCALE(k)  for  k  =  1,  ....  NPAR} 

This  convergence  test  is  roughly  equivalent  to  the  test  based  on 
the  maximum  relative  change  in  each  parameter  as  measured  by 
max{ l3£+i(k)-6£(k)|/l3£(k)|  fork  =  1,  ...,  NPAR} .  STOPP  is  not  a 
scale-dependent  value;  if  its  value  is  10-t+,  then  this  criteria 
will  he  met  when  the  first  four  digits  of  each  parameter  are  the 
same  at  two  successive  iterations,  regardless  of  the  size  of  the 
parameter  values. 

The  default  value  is  approximately  10~DIGITS/2 ,  where  DIGITS  is 
the  number  of  decimal  digits  carried  by  the  user's  computer  for  a 
single  precision  value  when  the  single  precision  version  of 
STARPAC  is  being  used,  and  is  the  number  carried  for  a  double 
precision  value  otherwise.  When  the  user-supplied  value  for  STOPP 
is  outside  the  interval  [0.0,  1.0],  or  when  STOPP  is  not  an 
argument  of  the  subroutine  CALL  statement,  the  default  value  will 
be  used. 

STOPSS  — >  The  stopping  value  for  the  convergence  test  based  on  the  ratio  of 
the  forecasted  change  in  the  residual  sum  of  squares,  As2£+l»  to 
the  current  residual  sum  of  squares,  s2 %  =  RSSCB^) /(NNZW-NPARE) . 
The  convergence  criterion  is  satisfied  if  certain  conditions  are 
met  [see  Dennis  et  al. ,  1981a]  and 

As2£+1/s2£  <  STOPSS. 

This  convergence  test  is  roughly  equivalent  to  the  test  based  on 
the  relative  change  in  the  residual  standard  deviation  between  two 
iterations  %  and  £+1  as  measured  by  (s^  -  s^+i) /$ %•  STOPSS  is  not 
a  scale-dependent  value;  if  its  value  is  10~5,  this  criteria  will 
be  met  when  the  first  five  digits  of  the  residual  sum  of  squares 
are  the  same  at  two  successive  iterations  regardless  of  the  size 
of  the  residual  sum  of  squares. 

The  default  value  is  approximately  the  maximum  of  10  and 
10-2* DIGITS/3 ^  Where  DIGITS  is  the  number  of  decimal  digits 
carried  by  the  user's  computer  for  a  single  precision  value  when 
the  single  precision  version  of  STARPAC  is  being  used,  and  is  the 
number  carried  for  a  double  precision  value  otherwise.  When  the 
user-supplied  value  for  STOPSS  is  outside  the  interval 
[10"DIGiTS,  0.1],  or  when  STOPSS  is  not  an  argument  of  the 
subroutine  CALL  statement,  the  default  value  will  be  used. 
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STP     The  vector  of  dimension  at  least  NPAR  that  contains   the  relative 

step  sizes  used  to  approximate  the  derivative  matrix  numerically. 
It  is  input  to  the  estimation  subroutines  and  returned  from  the 
step  size  selection  subroutines.  The  procedure  used  to  select  the 
default  values  is  described  in  §E.2.  For  the  estimation 
subroutines,  when  STP  is  not  an  argument  of  the  subroutine  CALL 
statement  or  when  STP(l)  <  0,  the  default  values  will  be  used  for 
all  of  the  step  sizes,  and  when  STP(l)  >  0,  values  of  STP(k)  <  0, 
k  =  2,  ...,  NPAR,  will  be  interpreted  as  an  input  error. 

VCV  < —  The  matrix  of  dimension  at  least  NPARE  by  NPARE  that  contains  the 
variance-covariance  matrix  of  the  estimated  parameters, 
approximated  as  designated  by  argument  IVAPRX.  The  parameters 
which  are  held  fixed  [see  argument  IFIXED]  are  not  included  in  the 
variance-covariance  matrix. 

The  approximation  of  the  variance-covariance  matrix  is  based  on  a 
linearization  of  the  model  in  the  neighborhood  of  the  solution; 
the  validity  of  the  approximation  depends  on  the  nonlinearity  of 
the  model.  This  approximation  may  be  extremely  inaccurate  for  a 
problem  with  a  highly  nonlinear  model. 

WT  — >  The  vector  of  dimension  at  least  N  that  contains  the  weights. 
Negative  weights  are  not  allowed  and  the  number  of  nonzero  weights 
must  equal  or  exceed  the  number  of  parameters  being  estimated.  A 
zero  weight  eliminates  the  corresponding  observation  from  the 
analysis,  although  the  predicted  values  and  standard  deviations  of 
the  predicted  values  are  still  computed.   [See  Appendix  A.] 

XM  — >  The  matrix  of  dimension  at  least  N  by  M  whose  j  column  contains 
the  N  values  of  the  jtn  independent  variable,  j  =  1,  ...,  M. 

Y  — >  The  vector  of  dimension  at  least  N  that  contains  the  dependent 
variable. 


E .   Computational  Methods 

E .1   Nonlinear  Least  Squares  Estimation 

E.l.a  The  Nonlinear  Least  Squares  Algorithm 

The  nonlinear  least  squares  estimation  subroutines  use  the  NL2S0L 
software  package  written  by  Dennis  et  al.  [1981a  and  1981b] .  The  observations 
of  the  dependent  variable,  which  are  measured  with  error,  are  iteratively  fit 
to  a  nonlinear  model  by  minimizing  the  sums  of  squares  of  the  errors  as 
described  in  §A.  The  iterations  continue  until  the  convergence  criteria  based 
on  the  change  in  the  parameter  values  or  in  the  residual  sum  of  squares  are 
satisfied  [see  §D,  arguments  STOPP  and  STOPSS] ,  the  maximum  number  of 
iterations  (or  model  subroutine  calls)  is  reached  [see  §D,  argument  MIT],  or 
the  iterations  are  terminated  due  to  singularity  in  the  model  or  false 
convergence.  All  but  the  first  of  these  stopping  conditions  may  indicate 
computational  problems  and  will  produce  an  error  report.  [See  Donaldson  and 
Tryon,  1983,  §D.5.] 
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Singular  convergence  means  that  the  model  contains  too  many  parameters, 
at  least  near  the  solution,  while  false  convergence  can  indicate  that  either 
STOPSS  or  STOPP  is  set  too  small  for  the  accuracy  to  which  the  model  and  its 
derivatives  are  being  computed,  or  that  there  is  an  error  or  discontinuity  in 
the  derivative.  Users  should  examine  their  models  to  determine  and  correct 
the  underlying  cause  of  singular  or  false  convergence. 

Iterative  procedures  for  solving  nonlinear  least  squares  problems  are 
discussed  in  Dennis  and  Schnabel  [1983] ,  Draper  and  Smith  [1981] ,  and  Kennedy 
and  Gentle  [1980] .  The  specific  procedure  used  in  STARPAC  is  as  follows.  At 
iteration  £+1  the  values  of  the  parameter  vector  $£+i  are  given  by 

H+l    =  H   ~    (V-W'D*  +   Si+    AjJ_1-D£T-W'e)tT 
subject  to  the  restriction  that 

NPAR 

{  I    [(0£+l(k>  "  Mk))/SCALE(k)]2}l/2  <  6£> 
k=l 


where 
NPAR 


is  the  number  of  parameters. 

is   the  vector  of  NPAR  estimated   parameter   values   from   the  I 
i  teration. 


th 


is  the  N  by  NPAR  matrix  of  the  partial  derivatives, 


D£(i,k)  =  8fi(xi,B£)/36)l(k) 


for  i  =  1,  . . .,  N 
and  k  =  1,  . . . ,  NPAR. 


W 


is  an  N  by  N  diagonal  matrix  of  user-supplied  weights, 
W  =  diag{wt^,  i  =  1,  ...,  N},  when  a  weighted  analysis  is  performed, 
and  is  the  identity  matrix  otherwise. 

is  an  approximation  to  the  exact  term  S^*  in  the  matrix  of  second 
order  terms  (Hessian)  of  the  Taylor  series  expansion  of  the  residual 
sum  of  squares  function, 

N 
S**(j>k)  =  I       [e£i-wti.(82e£(i)/3BJi(j)8B£(k))], 
i=l 


e£ 
SCALE 


for  j  =  1,  ...,  NPAR,  and  k  =  1,  ...,  NPAR. 

is  the  vector  of  the  N  residuals  from  the  previous  iteration. 

is  the  vector  of  the  scale  value  (i.e.,  the  typical  size)  of  each  of 
the  NPAR  parameters.  They  are  either  selected  by  STARPAC  or  are 
supplied  by  the  user  via  subroutine  argument  SCALE  [see  §D] . 
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6«  is  the  adaptively  chosen  diameter  of  the  trust  region,  i.e.,  the 
region  in  which  the  local  approximation  to  the  user's  model  function 
is  reliable.  At  each  iteration  i,  5^  is  computed  based  on  information 
from  the  previous  iteration.  At  the  first  iteration,  the  initial 
value,  6q,  is  supplied  by  argument  DELTA  [see  §D] ,  which  can  be  used 
to  control  the  change  in  the  parameters  permitted  at  the  first 
i  teration. 

A£  is  an  NPAR  by  NPAR  diagonal  matrix,  A£  =  diagj X^/SCALECk) ,  k  =  1,  ..., 
NPAR}  ,  where  Ag,  is  chosen  to  approximate  the  smallest  non-negative 
number  such  that  the  restriction  given  above  on  the  size  of  the  change 
in  the  parameters  is  satisfied. 

The  second  order  term  S^  ,  which  is  expensive  and  difficult  to  compute 
accurately,  is  important  only  if  it  is  large  compared  to  the  term  D^^.y.n., 
that  is,  when  the  residuals  are  large  or  the  model  is  highly  nonlinear.  When 
S^  is  large  compared  to  D£T»W»D£,  algorithms  which  ignore  it,  such  as 
Levenberg-Marquardt  or  Gauss-Newton,  may  converge  slowly.  The  NL2S0L 
algorithm  used  by  STARPAC,  however,  adaptively  decides  when  inclusion  of  this 
terra  is  necessary  for  reliable  results  and  uses  an  inexpensive  approximation 
to  Sj£  in  those  cases. 

The  matrix,  D,  of  partial  derivatives  of  the  model  with  respect  to  each 
parameter  is  either  computed  analytically  using  a  user-supplied  subroutine 
[see  §D,  argument  NLSDRV],  or  is  numerically  approximated  using  forward 
difference  quotients  as  described  in  §E.2.  When  the  derivatives  are 
approximated  numerically,  the  least  squares  solution,  especially  the 
variance-covariance  matrix,  can  be  sensitive  to  the  step  sizes  used  for  the 
approximation.  The  user  may  want  to  use  STARPAC  subroutines  STPLS  or  STPLSC 
to  recompute  the  step  sizes  at  the  solution  provided  by  the  estimation 
subroutines  to  assure  that  the  step  sizes  which  were  used  are  still 
acceptable.  If  there  is  a  significant  change  in  the  step  size,  the  least 
squares  solution  should  be  recomputed  with  the  new  step  sizes  from  the  current 
point.  In  addition,  if  the  estimation  subroutine  has  convergence  problems, 
the  user  may  want  to  recompute  the  step  sizes  with  the  most  recent  parameter 
values  to  see  if  a  change  in  the  curvature  of  the  model,  which  will  be 
reflected  as  a  change  in  the  optimum  step  sizes,  is  causing  the  problem. 

Dennis  et  al.  [1981a]  provides  a  detailed  description  of  the  algorithm 
used  in  STARPAC.  STARPAC  also  includes  the  subroutines  NL2S0L,  NL2SN0,  and 
NL2ITR,  which  they  reference,  and  which  can  be  used  as  documented  by  them. 
[See  Dennis  et  al. ,  1981b.] 


E.l.b   Computed  Results  and  Printed  Output 

From  the  Nonlinear  Least  Squares  Regression  Subroutines 

The  argument  controlling  the  printed  output,  NPRT,  is  discussed  in  §D. 
The  output  from  the  nonlinear  least  squares  estimation  subroutines  consists  of 
five  sections,  several  of  which  include  tables  summarizing  the  results.  In 
the  following  descriptions,  the  actual  table  headings  are  given  by  the 
underlined,  uppercase  phrases.  Results  which  correspond  to  input  or  returned 
subroutine  CALL  statement  arguments  are  identified  by  the  argument  name  in 
upper  case  (not  underlined) . 
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Section  1  provides  a  summary  of  the  initial  estimates  and  control  values. 
It  lists  the  following  information. 

•  The  initial  values  of  the  parameters,  PAR,  and  whether  they  are  to  be 
held  fixed  or  not,  IFIXED. 

•  The  scale  values,  SCALE. 

•  Either  the  step  sizes  used  to  approximate  the  derivatives  numerically, 
or,  when  user-supplied  (analytic)  derivatives  are  used,  the  results  of 
the  checking  procedure;  and  the  control  values  used  in  these 
computations  as  applicable  [see  §E.2.b  and  §E.3.b]. 

•  The  number  of  observations,  N. 

•  The  number  of  observations  with  nonzero  weights,  NNZW. 
9   The  number  of  independent  variables,  M. 

•  The  maximum  number  of  iterations  allowed,  MIT. 

«  The  maximum  number  of  model  subroutine  calls  allowed. 

•  The  two  convergence  criteria,  STOPSS  and  STOPP. 

»  The  maximum  change  in  the  parameters  allowed  at  the  first  iteration, 
DELTA. 

fl>  The  residual  sum  of  squares  computed  using  the  starting  parameter 
values. 

»  The  residual  standard  deviation  computed  using  the  starting  parameter 
values,  RSD. 

Section  2,  lists  selected  information  about  each  iteration  and  includes 
the  reason  the  iterations  were  terminated.  The  information  provided  for  each 
iteration  includes  the  following. 

•  The  iteration  number. 

9  MODEL  CALLS:  the  total  number  of  times  since  execution  began  that  the 
user's  model  subroutine  has  been  called,  not  including  calls  required 
to  approximate  the  derivatives  numerically. 

9  RSD:  the  residual  standard  deviation  computed  using  the  parameter 
values  from  the  current  iteration. 

•  RSS:  the  residual  sum  of  squares  computed  using  the  parameter  values 
from  the  current  iteration. 
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©  REL  CHNG  RSS:  the  relative  change  in  the  residual  sum  of  squares 
caused  by  the  current  iteration. 

•  FORECASTED  REL  CHNG  RSS:  the  forecasted  relative  change  in  the 
residual  sura  of  squares  at  the  current  iteration,  and  whether  this 
value  was  checked  against  STOPSS  (CHKD  =  Y)  or  not  ( CHKD  =  N) . 

•  REL  CHNG  PAR:  the  maximum  scaled  relative  change  in  the  parameters  at 
the  current  iteration,  and  whether  this  value  was  checked  against 
STOPP  (CHKD  =  Y)  or  not  (CHKD  =  N) . 

•  CURRENT  PARAMETER  VALUES:  the  estimated  parameter  values  resulting 
from  the  current  iteration. 

Section  3   provides   the   following   information   for  each  observation, 
i  =  1,  . . . ,  N,  based  on  the  final  solution. 

•  ROW:   the  row  number  of  the  observations. 

•  PREDICTOR  VALUES:  the  values  for  up  to  the  first  three  columns  of  the 
independent  variable  matrix,  XM ,  not  including  the  first  column  if  it 
is  constant. 

•  DEPENDENT  VARIABLE:   the  values  of  the  dependent  variable,  Y. 

©  PREDICTED  VALUE:   the  estimated  predicted  values,  PV,  from  the  fit. 

•  STD  DEV  OF  PREP  VALUE:  the  standard  deviations  of  the  predicted 
values,  SDPV. 

m   RESIDUAL:   the  error  estimates,  RES. 

•  STD  RES:   the  standardized  residuals,  SDRES. 

•  WEIGHT:  the  user-supplied  weights,  WT,  printed  only  when  weighted 
analysis  is  performed. 

Section  4,  displays  the  following  plots  of  the  standardized  residuals. 

»  The  standardized  residuals  versus  row  numbers. 

•  The  standardized  residuals  versus  predicted  values. 

•  The  standardized  residuals  versus  the  standardized  residuals  lagged  by 
one  observation,  that  is,  SDRES(i)  versus  SDRES(i-l) ,  i  =  2  . . . ,  N. 

©  The  normal  probability  plot  of  the  standardized  residuals. 

Section  5,  summarizes  the  following  information  about  the  final  parameter 
estimated  and  their  variances. 
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•  The  variance-covariance  matrix,  VCV,  of  the  estimated  (unfixed) 
parameters,  and  the  corresponding  correlation  matrix, 

rjk  =  VCV(j,k)/(vCV(j,j)  VCV(k,k))1/2      for  j  =  1,  ...,  NPARE 

and  k  =  1,  . . . ,  NPARE. 

•  PARAMETER:    the  final  value  of  each  parameter,  PAR(k),   k  =  1,  ..., 

NPAR. 

•  SD  OF  PAR:   the  standard  deviation  of  each  estimated  parameter, 
(vCV(k,k))1/2     for  k  =  1,  ...,  NPAR. 

•  RATIO:  the  ratio  of  each  estimated  parameter  to  its  standard 
deviation, 

RATIOk  =  PAR(k)/(vCV(k,k))1/2     f or  k  =  1 ,  ...,  NPAR. 

•  APPROXIMATE  9  5  PERCENT  CONFIDENCE  LIMITS:  the  lower  and  upper  95 
percent  confidence  limits  for  each  parameter,  computed  using  the 
appropriate  value  of  the  Student's  t  distribution  with  NNZW-NPARE 
degrees  of  freedom. 

•  the  residual  sum  of  squares,  RSS(3). 

•  the  residual  standard  deviation  at  the  solution,  RSD. 

•  the  residual  degrees  of  freedom,  NNZW-NPARE. 

•  an  approximation  to  the  condition  number  of  the  derivative  matrix,  D 
(the  Jacobian) ,  under  the  assumption  that  the  absolute  error  in  each 
column  of  D  is  roughly  equal.  The  approximation  will  be  meaningless 
if  this  assumption  is  not  valid;  otherwise  it  will  usually 
underestimate  the  actual  condition  number  by  a  factor  of  from  2  to  10 
[see  Dongarra  et  al. ,  1979,  p.  9.5]. 


NOTE:  The  standard  deviation  of  the  predicted  values,  the  standardized 
residuals,  the  variance-covariance  matrix,  the  standard  deviations  of  the 
parameters,  and  the  95  percent  confidence  limits  on  the  parameters  are  all 
based  on  a  linear  approximation  to  the  model  in  a  neighborhood  of  the 
solution;  the  validity  of  this  approximation  depends  on  the  nonlinearity  of 
the  model.  The  statistics  based  on  this  approximation  may  be  extremely 
inaccurate  for  a  problem  with  a  highly  nonlinear  model. 


E  .2   Step  Size  Selection 

E.2.a   The  Step  Size  Selection  Algorithm 

The  STARPAC  step  size  selection  subroutines  use  an  algorithm  developed  by 
Schnabel  [1982]  to  compute  optimum  step  sizes  for  approximating  the  partial 
derivatives  of  the  model  with  respect  to  each  parameter.    Briefly,   the 
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relative  step  sizes  selected  by  these  subroutines  are  those  which  produce 
forward  difference  quotient  approximations  to  the  derivative,  Df^j  that  agree 
reasonably  well  with  the  central  difference  quotient  approximations,  DC(j.  The 
central  difference  quotient  approximations  are  twice  as  accurate  but  also 
twice  as  expensive  to  compute.  Since  the  additional  accuracy  is  not  usually 
needed,  central  difference  quotient  approximations  are  not  used  by  the 
estimation  subroutines. 

The  number  of  reliable  digits  in  these  derivatives  is  a  function  of  the 
step  sizes  used  to  compute  them.  Given  properly  chosen  step  sizes,  the  number 
of  reliable  digits  in  Df^  and  Dccj  will  be  approximately  n/2  and  r\ , 
respectively,  where  n  is  the  number  of  reliable  digits  in  the  predicted 
values,  PV,  from  the  user's  model  subroutine.  For  example,  if  the  predicted 
values  are  computed  using  an  iterative  procedure  (such  as  quadrature  or  a 
solution  of  partial  differential  equations)  which  is  expected  to  provide  five 
good  digits,  then  n  would  be  five;  if  the  predicted  values  are  calculated  from 
a  simple  algebraic  expression  translated  directly  into  Fortran  code,  then  n 
would  (usually)  be  the  number  of  decimal  digits  carried  by  the  user's  computer 
for  the  results. 

The  relative  step  size  for  B(k),  k  =  1,  ...,  NPAR,  is  initially 

STP(k)  =  2(1CTNETA/Y)1/2     for  k  =  1,  ...,  NPAR, 

where 

Y       is   the  average  curvature   (estimated  by   STARPAC)  of   the  model  with 
respect  to  3(k) . 

NETA    is  the  argument  used  to  specify  the  number  of  reliable  digits,   n ,  in 
the  model  results. 

The  forward  difference  quotient  approximations  with  respect  to  3(k), 
k  =  1 ,  . . . ,  NPAR  are  then 

fi(xi,Bk)  -  £±(x±i$) 
Dfd(i,k)  =  for  1-1,  ...,  N, 

STP(k) -SCALE(k) •SIGN( B(k) ) 

where 

N       is  the  number  of  observations. 

f^      is  the  function  which  models  the  itn  observation. 


Xj  is  the  vector  of  the  values  of  the  M  independent  variables  at   the  i 

observation. 

B       is  the  vector  of  the  NPAR  parameter  values. 


th 
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8       is  a  vector  which  has   the   same   values   as   8,   except   that   the 
k1-   parameter  is  equal  to 

8(k)  +  (STP(k).SCALE(k).SIGN(3(k))). 

SCALE   is  the  vector  of  the  scale  value   (i.e.,  of  the  typical  size)   of  each 
of  the  parameters. 

SIGN    is  a  function  which  returns  the  sign  of  its  argument. 

The   central   difference  approximations   to  the  model  derivative  with 
respect  to  B(k),  k  =  1,  ...,  NPAR,  are 

fi(x±iB+h    -   f±(xlfB-k) 
Dcd(i,k)  =  for  1  =  1,  ...,  N, 

31/3.10"NETA/3.SCALE(k).SIGN(B(k)) 


where 

8       is  a   vector  which  has  the   same   values  as   8,  except   that   the  k*- 
parameter  is  equal  to 

8(k)  +  (31/3.10"NETA/3.SCALE(k).SIGN(8(k))). 

— k  t"  h 

8       is  a   vector  which  has  the   same   values  as   8,  except   that   the  kLU 
parameter  is  equal  to 

B(k)  -  (31/3.10"NETA/3.SCALE(k).SIGN(8(k))). 

The  relative  step  size  is  considered  acceptable  if,  for  at  least  N-a 
observations, 

|Dfd(i,k)  -  Dcd(i,k)|  <  min{lO-NETA/4)10-2|     for  i  =  1,  .....  N, 

where  a   is  the  number  of   observations   exempted  from  meeting   the  above 

acceptance  criterion  [see   §D,  argument  EXMPT] .    If   the   step  size   is  not 

acceptable,  it   is  adjusted  by  factors  of   10  until  the  condition  is  met,  or 

until  no  further  decrease  in  the  number  of  failures  can  be  made,   although  in 

no  case  will  the  selected  relative  step  size  be  greater  than  1.0  or  less  than 
j^-NETA^ 

Note  that  the  step  size  selection  subroutines  will  return  the  selected 
step  sizes  even  when  the  number  of  failures  exceeds  the  allowed  value;  this 
condition  will  be  noted  by  the  value  of  IERR.  The  detailed  printed  output 
should  always  be  examined  for  problems  discovered  by  the  step  size  selection 
subroutines. 
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E.2.b   Computed  Results  and  Printed  Output 

From  the  Derivative  Step  Size  Selection  Subroutines 

The  argument  controlling  the  printed  output,  NPRT,  is  discussed  in  §D. 

The  output  from  the  step  size  selection  subroutines  consists  of  a  summary 
of  the  input  and  control  values  and,  for  each  parameter,  the  selected  relative 
step  size,  the  number  of  observations  at  which  this  step  size  failed  the  step 
size  selection  criteria,  and  the   row  numbers  at  which  the   failures  occurred. 


E  .3   Derivative  Checking 

E.3.a   The  Derivative  Checking  Algorithm 

The  STARPAC  derivative  checking  subroutines  use  an  algorithm  developed  by 
Schnabel  [1982]  to  determine  the  validity  of  the  user-supplied  derivative 
subroutine.  The  user-supplied  derivative  subroutine  is  considered  correct  for 
a  given  row  i,  i  =  1,  . . . ,  N,  and  coefficient  B(k) ,  k  =  1,  ...,  NPAR,  if 


Dfd(i,k)  -  D(i,k)|  <  10-t  |D(i,k)| 


where 


'fd 


is  the  derivative  computed  by  the  user's  subroutine. 

is  the   forward  difference  quotient  approximation  to   the  derivative. 
[See  §E.2.a.] 

t  is  the  agreement  tolerance,  i.e.,  number  of  digits  of  agreement 
required  between  D  and  Df^,  which  must  be  less  than  or  equal  to  the 
number  of  good  digits  in  Df^.   [See  §D,  argument  NTAU.  ] 

When  the  agreement  tolerance  is  not  satisfied,  the  checking  subroutine 
attempts  to  determine  whether  the  disagreement  is  due  to  an  error  in  the 
user's  code  or  is  due  to  the  inaccuracy  of  the  difference  quotient 
approximation,  caused  either  by  high  curvature  in  the  user's  model  or  by 
significant  roundoff  error. 

The  derivative  checking  subroutines  each  check  only  one  row  of  the 
derivative  matrix.  The  user  should  examine  the  row  at  which  the  derivatives 
were  checked  to  ensure  that  some  relation  between  the  coefficients  and 
independent  variables,  such  as  a  zero-valued  independent  variable  or  a  factor 
(x^  -  3(k))  when  x^  =  $(k),  is  not  hiding  the  effect  of  an  incorrectly 
computed  derivative.  Checking  only  one  row  is  appropriate  since  the  same  code 
is  frequently  used  to  compute  the  model  function  and  derivatives  at  each  row 
i  =1,  . . . ,  N,  as  is  the  case  in  the  examples  shown  in  §F.  If  the  code  used 
to  express  the  model  function  and  derivatives  is  not  the  same  for  each  row, 
then  each  distinct  section  of  the  code  should  be  checked  by  making  multiple 
calls  to  DCKLSC  with  NROW  set  to  a  row  within  each  section.  [See  §D,  argument 
NROW.] 
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E.3.b   Computed  Results  and  Printed  Output 

From  the  Derivative  Checking  Subroutines 

The  argument  controlling  the  printed  output,  NPRT,  is  discussed  in  §D. 

The  output  for  the  derivative  checking  subroutines  consists  of  a  summary 
of  the  input  and  control  values,  and  the  results  of  the  derivative  checking 
test  with  respect  to  each  of  the  model  parameters,  3(k),  k  =  1,  •••,  NPAR. 
The  possible  test  results  are: 


OK  - 


•  The  user-supplied  derivative  and  the  numerical  derivative  agree  to  the 
required  number  of  digits. 


QUESTIONABLE  - 


•  The  user-supplied  derivative  and  the  approximated  derivative  agree  to 
the  required  number  of  digits  but  both  are  equal  to  zero.  The  user 
should  recheck  the  derivative  at  another  row. 

•  The  user-supplied  derivative  and  the  approximated  derivative  do  not 
agree  to  the  required  number  of  digits,  but  the  user-supplied 
derivative  is  identically  zero  and  the  approximated  derivative  is 
approximately  zero.  The  user  should  recheck  the  derivative  at  another 
row . 

•  The  user-supplied  derivative  and  the  approximated  derivative  disagree, 
but  the  user-supplied  derivative  is  identically  zero.  The  user  should 
recheck  the  derivative  at  another  row. 

o  The  user-supplied  derivative  and  the  approximated  derivative  disagree, 
but  the  validity  of  the  approximated  derivative  is  questionable 
because  either  the  ratio  of  the  relative  curvature  of  the  model  to  the 
slope  of  the  model  is  too  high,  or  SCALE(k)  is  wrong. 

•  The  user-supplied  derivative  and  the  approximated  derivative  disagree, 
but  the  validity  of  the  estimated  derivative  is  questionable  because 
the  ratio  of  the  relative  curvature  of  the  model  to  the  slope  of  the 
model  is  too  high. 


INCORRECT  - 


•  The  user-supplied  derivative  and  the  approximated  derivative  disagree, 
and  there  is  no  reason  to  question  the  accuracy  of  the  approximated 
derivative. 
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F .   Examples 

The  sample  programs  of  this  section  use  the  model  and  data  given  in 
example  one,  pages  428  to  441  of  Daniel  and  Wood  [1980];  the  model  is 

f±(xlt3)  =  6(l)-x1(l)3(2),     for  i  =  1,  ...,  N. 

In  the  sample  program  of  figure  F-la,  NLS  is  used  to  compute  the  least 
squares  solution  using  numerically  approximated  derivatives;  figures  F-lb 
through  F-lf  show  the  output  from  NLS. 

In  sample  program  of  figures  F-2a,  NLSD  is  used  to  compute  the  least 
squares  solution  given  analytic  derivatives;  figures  F-2b  through  F-2f  show 
the  output  from  NLSD. 

In  the  example  program  of  figure  F-3a,  STPLS  is  used  to  compute  the 
optimum  step  sizes  for  numerically  approximating  the  derivatives  with  respect 
to  each  of  the  parameters,  B(k) ,  k  =  1,  2.  Figure  F-3b  shows  the  output  from 
STPLS. 

In  the  example  program  of  figure  F-4a,  DCKLS  is  used  to  check  the 
validity  of  a  user-supplied  derivative  subroutine.  In  figure  F-4a,  the 
derivative  subroutine  has  been  intentionally  coded  incorrectly  in  order  to 
display  the  report  obtained  when  the  derivative  checking  subroutine  determines 
the  derivatives  are  incorrect,  and  the  starting  parameter  values  have  been 
chosen  in  order  to  display  the  report  obtained  when  the  test  results  are 
questionable.   Figure  F-4b  shows  the  output  from  DCKLS. 
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MAIN  PROGRAM! 


MODEL  SUBROUTINE! 


DATAl 


Nonlinear  Least  Squares  Estimation 
With  Numerically  Approximated  Derivatives 


PROGRAM  EXAMPL  (INPUT,  OUTPUT»  TAPES  ■  INPUT.  TAPE6  ■  OUTPUT) 


REAL  YUO),  XM(10»5),  PAR«5),  RES(IO) 
DOUBLE  PRECISION  DSTAK(2O0) 

EXTERNAL  NLSMDL 

COMMON  /CSTAK/  DSTAK 

OEFINE  VARIABLES  FOR  VARIOUS  DIMENSIONS 

LOSTAK  ■  200 
IXM  ■  10 

READ  NUMBER  OF  PARAMETERS 

STARTING  VALUES  FOR  PARAMETERS 

NUMBER  OF  OBSERVATIONS  AND  NUMBER  OF  INDEPENDENT  VARIABLES 

INDEPENDENT  AND  DEPENDENT  VARIABLES 

READ  (5,  100)  NPAR 

REAO  (5,  101)  (PARII),  I»1,NPAR) 

REAO  (5,  100)  N,  M 

READ  (5,  101)  MXH(I,J),  I-1.N)#  J"1,M),  (Yd),  I-1,N) 

PRINT  TITLE  AND  CALL  NLS  TO  PERFORM  NONLINEAR  REGRESSION 
WITH  NUMERICALLY  APPROXIMATED  DERIVATIVES 

WRITE  (6,  102) 

CALL  NLS  (Y,  XM,  N,  H,  IXM,  NLSMDL,  PAR,  NPAR,  RES,  LDSTAK) 

STOP 

FORMAT  STATEMENTS 

100  FORMAT  (215) 

101  FORMAT  (6F6.3) 

102  FORMAT  (1H1,  44HTHIS  IS  AN  EXAMPLE  OF  STARPAC  SUBROUTINE  NLS/) 
END 

SUBROUTINE  NLSMDL  (PAR,  NPAR,  XH,  N,  M,  IXM,  PV) 

SUBROUTINE  TO  COMPUTE  PREDICTED  VALUES  OF  OEPENOENT  VARIABLE 

REAL  PAR(NPAR),  XH(IXH,M),  PV(N) 

DO  10  1-1, N 

PV(I)  -  PAR(l)  *  XM(I,  1)  ♦•  PAR(2) 
10  CONTINUE 

RETURN 
END 


0.725  4.000 

6  1 

1.309    1.471    1.490    1.565    1.611    1.660 
2.138    3.421    3.597    4.340   4.662    5.660 


Figure  F-la 
Example  program  and  data  using  NLS 
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Nonlinear  Least  Squares  Estimation 
With  User-Supplied  Analytic  Derivatives 


MAIN    PR06RAMJ 


MODEL    SUBROUTINEi 


DERIVATIVE    SUBROUTINE 


OATAl 


PR06RAM  EXANPL  (INPUT,  OUTPUT,  TAPES  •  INPUT,  TAPE6  •  OUTPUT) 

REAL  Y(10)»  XM(10,9),  PAR(5),  REStlO) 
DOUBLE  PRECISION  DSTAM200) 

EXTERNAL  NLSMDL,  NLSDRV 

COMMON  /CSTAK/  DSTAK 

DEFINE  VARIABLES  FOR  VARIOUS  DIMENSIONS 

LDSTAK  •  200 
IXM  -  10 

READ  NUMBER  OF  PARAMETERS 

STARTING  VALUES  FOR  PARAMETERS 

NUMBER  OF  OBSERVATIONS  AND  NUMBER  OF  INOEPENOENT  VARIABLES 

INDEPENDENT  AND  DEPENDENT  VARIABLES 

READ  (9,  100)  NPAR 

READ  (9,  101)  (PAR(I),  I-1,NPAR) 

READ  (9,  100)  N,  M 

READ  (5,  101)  ((XM(I,J),  I-1.N),  J-1,M),  (Yd),  I-1,N) 

PRINT  TITLE  ANO  CALL  NLSD  TO  PERFORM  NONLINEAR  REGRESSION 
WITH  USER-SUPPLIED  DERIVATIVES 

WRITE  (6,  102) 

CALL  NLSD  <Y,  XM,  N,  M,  IXM,  NLSMOL,  NLSDRV,  PAR,  NPAR,  RES, 
*    LDSTAK) 

STOP 

FORMAT  STATEMENTS 

100  FORMAT  (215) 

101  FORMAT  (6F6.3) 

102  FORMAT  (1H1,  44HTHIS  IS  AN  EXAMPLE  OF  STARPAC  SUBROUTINE  NLSD/) 
END 

SUBROUTINE  NLSMDL  (PAR,  NPAR,  XM,  N,  M,  IXM,  PV ) 

SUBROUTINE  TO  COMPUTE  PREDICTED  VALUES  OF  DEPENDENT  VARIABLE 

REAL  PAR(NPAR),  XM(IXM,M),  PV(N) 

DO  10  I"1,N 

PV(I>  •  PAR(l)  ♦  XM(I,  1)  ♦•  PAR(2) 

10  CONTINUE 

RETURN 
END 

SUBROUTINE  NLSDRV  (PAR,  NPAR,  XM,  N,  M,  IXM,  D) 

SUBROUTINE  TO  COMPUTE  THE  PARTIAL  DERIVATIVE  (JACOBIAN)  MATRIX 

REAL  PAR(NPAR),  XM(IXM,M),  0(N,NPAR) 

DO  10  I-1,N 

0(1,1)  -  XM(I,1)  **  PAR(2) 

D(I,2)  ■  PAR(l)  *  XM(I,1)  **  PAR(2)  *  AL0G(XM(I,1) ) 
10  CONTINUE 

RETURN 
END 


0.729  4.000 

6  1 

1.309  1.471    1.490    1.969    1.611    1.680 

2.136  3.421    3.997    4.340   4.682    9.660 


Figure  F-2a 
Example  program  and  data  using  NLSD 
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Derivative  Step  Size  Selection 


MAIN  PROGRAM         PROGRAM  EXANPL  (INPUT,  OUTPUT*  TAPES  •  INPUT,  TAPE6  •  OUTPUT! 


REAL  *H!lO»5l,  PAB(5)*  STPI3) 
DOUBLE  PRECISION  DSTAM200) 

EXTERNAL  NLSNDL 

COMMON  /CSTAK/  DSTAK 


C 

C      DEFINE  VARIABLES  FOR  VARIOUS  DIMENSIONS 

C 

LDSTAK  -  200 

IXM  ■  10 
C 

C      READ  NUMBER  OF  PARAMETERS 
C  STARTING  VALUES  FOR  PARAMETERS 

C  NUMBER  OF  OBSERVATIONS  AND  NUMBER  OF  INDEPENDENT  VARIABLES 

C  INDEPENDENT  VARIABLES 

C 

READ  (5,  100)  NPAR 

READ  (5,  101)  (PAR(I),  1*1, NPAR) 

READ  (3,  100)  Ni  M 

READ  (3,  101)  ((XM(I,J>,  I-l,N),  J-1,H) 
C 

C      PRINT  TITLE  AND  CALL  STPLS  TO  SELECT  STEP  SIZES  FOR 
C      APPROXIMATING  DERIVATIVES 
C 

WRITE  (6,  102) 

CALL  STPLS  (XN,  N,  H»  IXM,  NLSMDL,  PAR,  NPAR,  LDSTAK,  STP) 
C 

STOP 
C 

C      FORMAT  STATEMENTS 
C 

100  FORMAT  (215) 

101  FORMAT  (6F6.3) 

102  FORMAT  (1H1,  46HTHIS  IS  AN  EXAMPLE  OF  STARPAC  SUBROUTINE  STPLS/) 
END 

MODEL  SUBROUTINE!         SUBROUTINE  NLSMDL  (PAR,  NPAR,  XM,  N,  M,  IXM,  PV) 
C 

C      SUBROUTINE  TO  COMPUTE  PREDICTED  VALUES  OF  DEPENDENT  VARIABLE 
C 

REAL  PAR(NPAR),  XM(IXN,M),  PV(N) 
C 

DO  10  I  ■  1,  N 

PV(I)  •  PAR(l)  •  XM(I,  1)  •*  PAR(2) 
10  CONTINUE 
C 

RETURN 
END 

DATAl       2 

0.725  4.000 

6  1 

1.309    1.471    1.490    1.565    1.611    1.680 


Figure  F-3a 
Example  program  and  data  using  STPLS 
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Derivative  Checking 


MAIN  PROGRAM! 


MODEL  SUBROUTINE! 


DERIVATIVE  SUBROUTINE! 


PROGRAM  EXAMPL  (INPUT*  OUTPUT,  TAPES  ■  INPUT,  TAPE6  -  OUTPUT) 

REAL  XM(10,5),  PAR(5) 
DOUBLE  PRECISION  DSTAK(ZOO) 

EXTERNAL  NLSMOL,  NLSORV 

COMMON  /CSTAK/  DSTAK 

DEFINE  VARIABLES  FOR  VARIOUS  DIMENSIONS 

LDSTAK  ■  200 
IXM  «  10 

READ  NUMBER  OF  PARAMETERS 

STARTINS  VALUES  FOR  PARAMETERS 

NUMBER  OF  OBSERVATIONS  ANO  NUMBER  OF  INDEPENDENT  VARIABLES 

INDEPENDENT  VARIABLES 

READ  (5,  100)  NPAR 

READ  (3,  101)  (PAR(I),  I- 1, NPAR) 

READ  (5,  100)  N,  M 

READ  (5,  101)  ((XM(I,J),  I-1,N),  j>l,M) 

PRINT  TITLE  AND  CALL  DCKLS  TO  PERFORM  DERIVATIVE  CHECKING 

WRITE  (6,  102) 

CALL  DCKLS  (XM,  N,  M,  IXM,  NLSMDL,  NLSORV,  PAR,  NPAR,  LDSTAK) 

STOP 

FORMAT  STATEMENTS 

100  FORMAT  (215) 

101  FORMAT  (6F6.3) 

102  FORMAT  (1H1,  46HTHIS  IS  AN  EXAMPLE  OF  STARPAC  SUBROUTINE  DCKLS/) 
END 

SUBROUTINE  NLSMDL  (PAR,  NPAR,  XM,  N,  M,  IXM,  »V) 

SUBROUTINE  TO  COMPUTE  PREDICTED  VALUES  OF  DEPENDENT  VARIABLE 

REAL  PAR(NPAR),  XM(IXM,M),  PV(N) 

DO  10  I-1,N 

PV(I)  "  PAR(l)  ♦  XM(I,1)  ♦♦  PAR(2) 

10  CONTINUE 

RETURN 
END 

SUBROUTINE  NLSORV  (PAR,  NPAR,  XM,  N,  M,  IXM,  0) 

SUBROUTINE  TO  COMPUTE  THE  PARTIAL  DERIVATIVE  (JACOBIAN)  MATRIX 

DERIVATIVE  WITH  RESPECT  TO  FIRST  PARAMETER  HAS  BEEN  CODED 
INCORRECTLY  (WITH  MULTIPLICATION  RATHER  THAN  EXPONENTIATION) 
FOR  DEMONSTRATION  PURPOSES 

REAL  PAR(NPAR),  XM(IXM,M),  DIN, NPAR) 

DO  10  I-1,N 

D(I,1)  ■  XM(I,1)  •  PAR(2) 

0(1,2)  •  PAR(l)  »  XH(I,1)  ♦*  PAR(l)  *  ALOG(XN(I,D) 
10  CONTINUE 

RETURN 
END 


DATA! 


0.000  4.000 

6     1 
1.309  1.471  1.490  1.365  1.611  1.680 


Figure  F-4a 
Example  program  and  data  using  DCKLS 
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Appendix  A 
WEIGHTED  LEAST  SQUARES 


Weighted  least  squares  can  be  used  to  eliminate  observations  from  the 
analysis,  and  to  compensate  for  unequal  variances  in  the  observational 
errors. 

Observations  can  be  eliminated  from  the  analysis  by  using  weight  values 
consisting  only  of  zeros  and  ones.  This  will  produce  the  same  results  as 
performing  an  unweighted  analysis  with  the  zero-weighted  values  removed  except 
that  the  predicted  values,  the  standard  deviations  of  the  predicted  values, 
and  the  residuals  of  the  zero-weighted  data  are  computed.  There  are  two  main 
reasons  for  weighting  observations  zero.  The  first  is  to  obtain  the  predicted 
values  and  their  standard  deviations  for  a  set  of  independent  variables  not 
included  in  the  observed  data.  (This  is  done  by  assigning  any  arbitrary  value 
to  the  dependent  variable  of  the  desired  set  of  independent  variables,  and 
then  weighting  these  values  zero.)  The  second  reason  is  to  allow  easy 
examination  of  the  effect  of  outliers  and  influential  data  points.  Outliers 
often  appear  as  large  values  in  residual  plots.  Careful  checking  of  the  data 
often  leads  to  confirmation  that  the  data  are  in  error,  and  sometimes  to  a 
correction.  When  a  cause  for  suspicious  data  cannot  be  found,  it  may  be 
advisable  to  compare  the  analysis  with  and  without  the  questionable  data. 
Caution  is  in  order  if  the  estimates  or  conclusions  are  highly  sensitive  to  a 
small  amount  of  suspicious  data.  Data  that  have  a  very  high  influence  on  a 
fitted  curve  may  not  result  in  large  residuals,  however,  even  if  they  are  in 
error.  In  fact,  extremely  influential  observations  may  force  the  fitted  curve 
to  be  very  close,  leading  to  very  small  residuals.  It  is  therefore  desirable 
to  identify  influential  observations  and  to  compare  the  results  obtained  with 
and  without  these  points.  Several  methods  for  detecting  influential 
observations  are  discussed  in  Bement  and  Williams  [1969],  Cook  [1977],  Hoaglin 
and  Welsch  [1978],  and  Belsley  et  al.  [1980]. 

Using  weights  to  compensate  for  unequal  observational  error  variances  is 
not  as  straightforward  as  using  zero  weights  to  eliminate  observations  from 
the  analysis.  When  the  variances  of  the  observational  errors,  e^,  are  not 
equal,  the  unweighted  least  squares  estimates  remain  unbiased  but  do  not  have 
minimum  variance.  Minimum  variance  estimates  are  obtained  by  using  weights 
wt.j  =  1/Variance[ e^]  when  the  error  variances  are  known .  If  weights  must  be 
estimated,  they  should  be  based  on  at  least  10  degrees  of  freedom  [see  Bement 
and  Williams,  1969].  In  practice,  however,  weights  are  derived  from  theory, 
or  obtained  from  the  data  being  fit,  and  either  of  these  methods  can  do  more 
harm  than  good.  When  the  need  for  weights  is  suspected  and  the  error 
variances  are  not  known,  first  fit  the  data  using  unweighted  least  squares; 
analysis  of  the  residuals  may  confirm  the  need  for  weighting  and  may  also 
provide  estimates  for  the  weights  themselves.  If  the  need  for  weights  is 
confirmed,  then  a  statistician  should  be  consulted  to  assist  in  selecting  the 
weights  and  in  interpreting  the  results. 
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Appendix  B 

ESTIMATING  THE  NUMBER  OF  RELIABLE  DIGITS 
IN  THE  RESULTS  OF  A  FUNCTION 


The  number  of   reliable  digits,   n,   in  the  results  of  a  real  valued 
function,  g((3),  can  be  estimated  in  most  cases  by  evaluating 

n  =  -log10(    max    (lg(Bj)  -  [a  +  j»b] I}] 
J — 2 2        fiOOl 

where 


6 J      is  the  vector  of  the  NPAR  parameters  of  the  function  given  by, 

BJ(k)  =  B(k)  +  j.lO-(DIGITS/2).B(k)     for  j  =  -2,  ...,  2, 

and  k  =  1,  . . .,  NPAR, 

where 

DIGITS  is  the  number  of  decimal  digits  carried  by  the  user's  computer 
for  a  single  precision  value  when  the  single  precision  version  of 
STARPAC  is  being  used  and  is  the  number  carried  for  a  double  precision 
value  otherwise. 

2 
a       =  (0.20).  y    g(3J). 
J— 2 

2 

b  =  (o.io).  y     j.  g(0j). 

j— 2 

This  procedure  may  underestimate  the  number  of  reliable  digits  if  g(3)  is 
extremely  nonlinear.  A  more  elaborate  and  more  robust  procedure  is  described 
in  Gill  et  al.  [1981]  . 
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